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Abstract 



We consider the statistical distribution of zeros of random meromorphic 
functions whose poles are independent random variables. It is demonstrated 
that correlation functions of these zeros can be computed analytically and 
explicit calculations are performed for the 2-point correlation function. This 
problem naturally appears in e.g. rank-one perturbation of an integrable 
Hamiltonian and, in particular, when a <5-function potential is added to an 
integrable billiard. 
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I. INTRODUCTION 



The investigation of statistical properties of quantum energy levels of a given system is a 
long-standing problem (see e.g. @). According to the accepted conjectures energy levels 
of integrable systems behave as independent random variables (i.e. they obey the Poisson 
statistics) [|j] and those of generic chaotic systems follow the random matrix predictions ||. 
The proof of these conjectures in the full generality is without doubt quite difficult and is 
still lacking though partial results (concerning mostly integrable models) are available (see 
e.g. || and references therein). 

But there are systems which are neither integrable nor completely chaotic for which 
quantum energy levels are defined by an equation 

f(E) = (1) 

with a well defined (and simple) function f(E). 
In the case of polynomial equation 

N 

f(E) = £ a n E n (2) 

n=0 

has been considered and statistical properties of solution of f(E) = have been calculated 
provided a n be independent random variables. 

The purpose of this paper is to consider the case of random meromorphic functions of 
the form 

N r- 

f(E) = P(E) + £ (3) 

j=l E ~ e 3 

where P(E) is a polynomial and ej, r 3 - are, correspondingly, poles and residues of f(E). 

The natural example leading to the quantization condition in this form is the perturbation 
of a Hamiltonian by rank-one perturbation. If Hffj) is an unperturbed Hamiltonian then the 
Hamiltonian after perturbation is 

= Hj$ + v*, (4) 

where t> M is a perturbation vector. 

Solutions of the 'Schrodinger' equation 

H^ u = E^ (5) 

can be expressed through solutions of unperturbed equation 

H^\n)=e n ^\n) (6) 

as follows 

^ = £<v/f(n), (7) 
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where (up to a factor) 



< u|V (0) (n) > ^ 



E - e n 

provided new eigenvalues, E, obey the following quantization condition 



< v\ip { - Ql \n) > 



2 



E E-e =L (9) 

n 1 J c r?. 



Here < v\^°\n) >= E^V/^H- 

This equation has the form of Eq. (|3|) with P(i?) =const while unperturbed energy levels 
play the role of poles, and the residues are projections of unperturbed wave functions in the 
directions of perturbation vector 

r n = \< v\^ Q) (n) > | 2 . (10) 

The addition of a 5-function potential 

V(x) = X5(x-xo) (11) 

corresponds exactly to a rank-one perturbation. In this case (see e.g. ||, H) Eq. @ takes 
the form 

AE MW = 1 , (12) 

where ip^(x) and e n are eigenfunctions and eigenvalues of the problem without the 5- 
function potential. 

Another model which leads to similar equations is the Bohr-Mottelson model JTOf which 
describes the interaction of one level (denoted below by index 0) with all other levels. The 
model is defined by the Hamiltonian 

H = H + V, (13) 

where the interaction potential has non-zero matrix elements only between the chosen level 
and all other levels 

V 0i = V i o, ^oo = ^ = 0. (14) 



The energy levels of the Hamiltonian (|13D obey the equation [10 



W 0j 



|2 



Ef^-(^-^o) = 0, (15) 

which is again of the form of Eq. @ with linear polynomial part. 

A quite natural question appears: What is the statistical distribution of new eigenvalues 
(i.e. solutions of Eq. (|3])) provided that statistical distributions of poles and residues are 
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known? In |TT| it was proved that, if unperturbed system is described by random matrix 
theory, the distribution of new eigenvalues will also be of random matrix type. 

The main purpose of this paper is to compute analytically the statistical distribution of 
solutions of Eq. (|3]) when the poles, ej, are independent random variables (i.e. they obey 
the Poisson statistics). We shall show that in this case the resulting statistics exhibits the 
level repulsion and differs from known distributions. 

The plan of the paper is the following. In Section |TT] the general formalism is described. 
In Section |T| the calculation of the mean density is presented. In Section [IV| the 2-point 
correlation function is computed when all residues, Vj, in Eq. (^) are the same. Generalization 
to different residues is disscused in Section |V11| . As the exact expression of the 2-point 
correlation function is cumbersome, in Section |V| the series expansion of the results is given. 
In Section |VT] the limiting behavior of the 2-point correlation function for small and large 
energy difference is obtained without the knowledge of the exact solution. The details of 
the calculation of a certain important integral are presented in Appendix. 



II. GENERAL FORMALISM 

We consider the most interesting case of Eq. (|3p when the mean separation of the poles 
is much smaller than a characteristic scale of polynomial P(E). Under such condition this 
polynomial can be considered as a constant and after dividing by it Eq. (HD takes the form 

j=l ^ e J 

Our goal is to find the statistical distribution of solutions, E, of this equation provided r 3 - are 
constants and N numbers Cj are independent random variables with a common distribution 
dfj,(e) which for simplicity we choose as follows 

da(e) = l^v deii - W ^ e > ^ W > . (17) 
1 otherwise 

As the density of these poles is a constant they can be considered as eigenvalues of a 2- 
dimensional integrable billiard and we shall often call them energy levels (or unperturbed 
energy levels). All our calculations remain also valid in more general case when the mean 
density of poles is not a constant but is not changed noticeable in the scale of the mean 
pole separation (e.g for 3-dimensional integrable models). The only difference is that N/2W 
below should be substituted by the local mean density of poles, p (see the end of Section 



In general, if one is interested in solutions of equation 

f(Xn) = 0, (18) 

it is often convenient to express the exact density of such solutions 

p{x)=Y,8(x-x n ), (19) 
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in the following manner 



p(x) = 5(f(x))\^\. (20) 



The main advantage of such representation is the possibility of calculating the statistical 
distribution of roots, x n , directly from statistical distribution of coefficients of f(x). E.g. 
this method has been used for deriving the distribution of roots of random polynomials ]7] . 
In our case 



v u " u 



Representing the 5-function as the Fourier integral (i.e. considering the characteristic func- 
tion of the roots) one gets 

^ - £ d £ ^ " 1}) £ (22) 

It is this representation of the exact density that we shall use throughout the paper. 

As all 6j are considered as independent random variables this expression can be rewritten 
in the form 

^) = /^-'«nexp(^)i:^, (23, 

where all factors are also independent random variables which clearly permits to find all 
mean values by straightforward integration. 



III. MEAN DENSITY 

Let us start with the calculation of the mean density 

< P( E ) >= / H / Me,) exp(^) £ (24) 



3=1 



The integrals can be transformed in the following way 



where 



and 



< P(E) >= I pe'™ r k g(r k a)) ]J f(r,a)) , (25) 

J Z7F \k=l i+k ) 



f(a) = j d/i(e)exp(i^-—), (26) 
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/I a <9 2 

Me) W^W eMl ^ ] = ~d^ f{a) - (27) 



Let us rewrite the expression for f(a) in the form 



f(a) = l-^I(a), (28) 



where 

rW 



1(a) =[ de(l-exp(i-^)), (29) 
J-w E — e 



As 



^ = W 1 ^' < 30 > 

it is necessary to compute only 1(a). 

Though the above steps are exact for finite N, the most interesting case is the case 
iV — > oo. In this limit only small values of a are important (a ~ and it is necessary 

to take into account in 1(a) only terms linear in a. 

Due to the singular character of the integral 1(a) (p9|) one cannot just expand the 
integrand in power of a. If E belongs to the support of the measure, — W < E < W, the 
change of variable 

*=^b (31) 

reduces the integral for 1(a) (|29|) to a sum of two integrals 

/ r oo f (E-W)- 1 \ df 

/(«)= / + / (32) 

which can be transform as follows 

^=\L-L w) -M- e ^- (33) 

The first integral equals ir\oe\ and in the second integral one can safely use perturbation 
theory in a. The final result is 

, • i W-E 2 W _ SN 

1(a) = ir\a\ + la In — - - + a 1 ^ _ ^ 2 + 0(a A ), (34) 



and 



9(a) = i«(a) + ^i^j. (35) 



For small values of a 
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e~ ia ft f(rja) = exp(-^v(n\a\ + (36) 



where t> plays the role of a 'bare' coupling constant, 



JV j=l 

and u' is a 'renormalized' coupling constant 

1 2W , W-E 

— = + ln . (38) 

v' Nv W + E K ' 

The necessity of renormalization for such type of equations is well known when a ^-function 
potential is added to a <i-dimensional system with d > 2 (see e.g. || and Eqs. (|98]and ( |131| )) 
where it is connected with one-parameter self-adjoint extension of a singular Hamiltonian. 
Physically the renormalization means that the limit of infinite small size impurity is not 
uniquely defined and depends on internal details of the scatterer. All physically measurable 
quantities (like the cross-section) depend only on renormalized coupling constant, v' . The 
bare coupling constant, v, is not observable and can be arranged to produce any v'. When a 
specific model of small-size scatterer is considered (e.g. a hard disk with a small radius) one 
gets a concrete form of the bare (and renormalized) coupling constant. Below we consider 
the most interesting case when renormalised coupling constant is assumed to be independent 
of N (or energy). All other limits can be derived from this one. Note that in our calculations 
the appearance of such renormalization (i.e. the fact that the bare coupling constant, v, 
and the renormalization factor log(W — E)/ (W + E) appear only in the combination (|38|)) 
is automatic. 

Finally, when —W<E<W the density of state is the sum of two terms 

N 2W 

Pin(E) = — ~ (W 2_ E 2 ) ^2 +l/v ,2y ( 39 ) 

As N is assumed to be large, the first term dominates and the mean density of levels is 

"=w- (40) 

as it should be. 

When E is beyond the interval [— W, W] the calculation is simpler as in this case there 
is no singularity on the contour of integration and one can simply expand the integrand of 
1(a) on series of a 

1{a) = ia ln §TW + a2 E^W> + 0{a3) - (41) 



Therefore 



^ E)=] M^--h l) - (42) 
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where 



, w = h £±E_Iv=«?Z_ (43) 

vy ! E-W v r (E 2 -W 2 )N V ; 

When N — ► oo, a — ► and 

p(E)^5(E-E c ), (44) 
where £7 C is a root of equation 0(£ , c ) = 

^Wcoth^. (45) 

These results correspond exactly to what one attends from the simple geometrical picture 
of roots of Eq. fllBT) . The poles, ej, divide the real axis into N +1 intervals. Due to the pole 
behavior each interval contains one of solutions, E. There is only one eigenvalue outside of 
the support of the initial measure and all other (N—l) eigenvalues are distributed practically 
uniformly inside the initial interval [— W, W\. The second term in Eq. (EDI) is a smooth bump 
which is necessary to insure that 



/w 
Pm {E)dE = N - 1, 
-w 



which can easily be checked by noting that 2W/(E 2 — W 2 ) = d(l/v')/dE. 

In Eq. (|TTD we have assumed the particular form of the distribution of dp(e) but the 
results will be valid for any form of this measure (provided that it is not changed noticeably 
in the scale of the mean distance between levels) with the substitution N/2W — > p, E+W — ► 
E — E min and W — E — » E max — E where p is the local mean density of unpertubated levels, 
E min and E max are minimal and maximal values of levels included in the sum ([16]) . 



IV. 2-POINT CORRELATION FUNCTION 

Using the previously discussed method one can compute higher correlation functions 
as well. We consider here the calculation of the 2-point correlation function, R2(E\, E2), 
defined in the standard way 

R2(E 1 ,E 2 )=<p(E 1 )p(E 2 )>, (46) 

where < . . . > denotes the mean value over all random variables. 

For clarity we consider first the case where all residues are equal, = v. This case 
appears e.g. when a 5-function potential is added to a rectangular billiard with periodic 
boundary conditions (see Eq. fl58|)). More general case with different will be considered 
shortly in Section |VII| . 



When all residues are the same our defining equation takes the form 



N 1 x 

S^ = «- (47) 
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and the 2-point correlation function can be expressed as follows 

, , f da x da 2 ,.^U f ai t a 2 
R 2 {E 1 ,E 2 ) = < / 2 exp(z^(- + - 



N i 

V - e -^(«i+a 2 ) > f48l 

J^ =1 - e ki nE 2 - e k2 y ■ { > 



After a simple algebra this expression can be transformed to 

da\da 2 



R 2 (E U E 2 ) = J -^-^[Nifiaua^f-'gia^a,) (49) 
+N(N - l)(/(ai, a 2 )) N ~ 2 i)i(ai, a 2 )^ 2 («i, a 2 )] exp(--(ai + oc 2 )), 



where 



/(ai,a 2 ) = / d/x(e) exp(i— — hi 

j Et — e 



E 1 -e E 2 -e 
g(a\,a 2 ) = I dfi(e)exp(i— — 



E 1 -e E 2 -e / (E 1 -e) 2 (E 2 -e] 



V>i(ai,a 2 )= I dfi(e)exp(i ai + i " 2 ) 1 (50) 

£/x — e E 2 — e (Ei — e) z 

. Oil ■ Oi 2 1 



d/x(e) exp(i— h i 

Ei — e 



Ei-e E 2 -e'(E 2 -e) 2 

We shall be interested in the distribution of eigenvalues inside the interval [-W, W] and 
therefore shall assume that both arguments E\ and E 2 belong to this interval. 
Let us as denote 

/(an, a 2 ) = 1 - 7^7 J ( a i> "2), (51) 

where 

f w oil ce 2 

I(ai,a 2 )= (l-exp(z— ))de. (52) 

J-w Ei — e E 2 — e 

Other functions are expressed through J(a 1 ,a 2 ) as follows 

^ (ai ' a2) = -2W^| /(ai ' a2) ' 



ZWdaj 
1 8 2 



^i(«i,a 2 ) = OTJAa 2 /(ai,a 2 ), (53) 



^ 2 (ai,a 2 ) = ^^| J ( ai ' a2 )- 
The integral (|52"D which defines J(ai,a 2 ) can be split into three terms 

/(«i, a a ) = f r - f ^ - (1 - exp(z-^- + z-^ 2 -))^. (54) 
\j—oo J-oo Joo l hii — e rLi 2 — e 



In the first integral (which we denoted by J(ai,a 2 )) singular points E\ and E 2 are on the 
contour of integration. In the second and the third ones there are no singularities and they 
can be computed in perturbation theory on ot\ and a 2 . In the later integrals we will see that 
one needs only linear in a terms and 

W - Ex W - E 2 

I(a x , a 2 ) = J(ai, a 2 ) + i(a x In W + E + «2 In w + ). (55) 

It is the calculation of the first term which is difficult. The details of this calculation are 
given in Appendix. The final result for J(ai,a 2 ) is the following 

Ot 2 Oil 

J(«i, a 2 ) = n(ai + a 2 )sgn(a 2 ) - 7r[z(ai + a 2 )G( , — ) (56) 

D U! 



+(a 2 Jo(0 + W-ocia 2 Ji(0) ex P(^ )](sgn(ai) - sgn(a 2 )), 

UJ 

where uj — E\ — E 2 , £ = ^\f—a\a 2 and 

roo 

G(x,y) = e ty J (2Vyi)e u dt. (57) 

J X 

The symmetry relations 

J(a 2 , «i) = J*(«i, a 2 ), 
J(—ai, — a 2 ) = J*(ai, a 2 ), (58) 
J(— a 2 , — at) = J{oii,a 2 ), 

are also useful. We are interested in the situation when the difference of energies, to = 
Ei — E 2 , is of the order of the mean distance between the levels 

U = n™, (59) 

and dimensionless frequency Q is a constant. In this case one can check that the important 
values of a will also be of the order of 1/N which explains why we have restricted the expan- 
sion only up to linear terms. Other simplification comes from the fact that in perturbation 
theory terms fl5"5]) one can put Ei = E 2 after which they only depend on the sum ai + a 2 . 
In the limit of large N one obtains (see Appendix) 

N ~ 

f (ai,a 2 ) = exp(-— I(a 1 ,a 2 )), (60) 

11.9 d ai — a 2 

9{ai > a2) = 2W-^da~i - ^ ^^^Woii, a 2 )}, (61) 

1 ai — a 2 d 

V>i{oci, a 2 ) = TyyTF ex P(« ^(ai, a 2 ), (62) 

2W u oaii 

1 Gii — 9 

ip 2 (ai,a 2 ) = -— exp(i — - — )— $(ati, a 2 ). (63) 
Here we introduce 
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I (aii, a 2 ) = J{ol\, 012) + («i + "2)^7, (64) 



where f ' is the renormalized coupling constant as in (RBI) and 



$(ai,a 2 ) = 27rJ (( — v/-aia 2 )6(-Q;ia2)sgn(ai). (65) 



Therefore 



R 2 (uj) = f daida2 \Nf N ~^(^- d -)\e^<$>] 

- N{N - Vf-v^A-wJLw. m 

It is convenient to integrate the first term by parts 

[ r-v^$(— - —)j = re?) 

N -l r ^ Ar _ 2e2i £i_i!2 $2 



2W 

Substituting this expression into the previous equation one obtains 

The second useful form can be derived by the following transformation of the second term 

= ^ + h& ' + >'«*• (69) 



where 



* = 2i^^. (70) 



CO 

Combining these two expressions one gets 

^2 (w) = — ( \ ' / daida 2 e-™ T — + -(- — $ 2 e*. 71 

It is easy to check that under the scale transformation (assuming A > 0) 

c<j — > \u, ati — > Aa«, (72) 
the pre-factor does not change and / — > XI. Therefore after the transformations 

n = (73) 
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and 

W = ^^r 2 (n), (74) 

plus the corresponding change of a the dependence of N will disappear and after the sub- 
stitution 

a.i = VLoii (75) 
the resulting expression for the 2-point correlation function takes the form 

r . m = _/^ e --,_|_ _ A )]42e , ( «-«, (76) 



where 



and from Eq. (|191|) 



J(ax,a 2 ) = J(a 1 ,a 2 ) + i(a x + a 2 )—— (77) 

ZlTV 



J{a 1 ,a 2 ) = -gfai + «2)sgn(a 2 ) 



+ -[(ax + a 2 )G(ax, -a 2 ) - (ia x J (2 v / -«i«2) (78) 
+ v /3 ^ 2 "Ji(2v /3 ^ 2 "))e i(ai - a2) ](sgn(a 1 ) - sgn(a 2 )). 
When Q — > it is convenient to perform the integration by part 

r 2 (0) = - / daxda 2 e- 2 ^ 3 [%^-^- + ^ e ^-^ e M^) (79) 
J 2 a ax oax 87H 

and take into account only the linear in Q term 

r 2 (Q) -> QA, (80) 

where 

A = y d ai (ia 2 $ 3 e 3i(ai - a2) , (81) 



As in the region aia:2 < $ = 27rJ (2v /— CKiaJ), after the change of variables 



ax 



e = 2v^^, 77= L -, (82) 

one gets 

A= .^r (4m r%w^ + c,.. (83) 

Z JO JO Tj 

The integral over rj equals mH^ (3£) (see (|169 )) and the final expression for A is 
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/■oo 

A = 3vr 2 lim/ ^o((3 + e)0^m- (84) 

We write here (3+e) (where e is proportional to Q) as this integral is a discontinuous integral 
and its value when e = is a half of the value for e — > 0. The last value can be computed 
using the following integral ( [[HJ, p. 414) 

/ II J {a n t)tdt = - < i ' / x ., -. ; (85) 

where is the full elliptic integral of the second kind 

x = A (86) 

and 

4 

16A 2 = n(«i + a 2 + a 3 + a 4 -2a n ). (87) 

n=l 

If the left hand side is negative the above integral equals zero. 
In our case A — ► and K(0) = tt/2, therefore 

/•OO I 

Mmjf o ^o((3 + eX)J 3 (0^=— y=. (88) 



Hence 



A = ^5 « 2.72.... (89) 



Note that the slope at the origin is independent on the coupling constant and differs from 
the prediction of the Gaussian Orthogonal Ensembles of random matrices (^(O) — > (tt 2 /6)Q 

0)- 

To find the asymptotics of the 2-point correlation function when Q — > oo it is convenient 
to use Eq. (pq). After rescaling of this expression one obtains (the constant term comes 
from the 5-function contribution of derivatives) 

{/•oo rO 
J dai J da 2 [Jl(2yJ-a 1 a 2 ) + J 1 2 (2 v /-aia 2 )] 

x exp (-2nnj + 2i(a x - a 2 )) + c.c] . (90) 

When Q — > oo the dominant contribution comes from region of small a. Taking into account 
that when a — > 

J^-( ai -a 2 )+E — , (91) 
one concludes that the corresponding asymptotics of the 2-point correlation function is 

r2(n > - 1 + n^li/v-y (92) 
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Note the absence of oscillation on large typical for standard random matrix ensembles. 

To check the above results we compute the statistical distribution of energy levels of a 
rectangular billiard with a 5-function potential inside (sometimes called the Seba billiard 

!)■ 

For a rectangle of sides a and h solutions of the Schrodinger equation 

(e n - A)Vft(£) = (93) 
in 2 dimensions with periodic boundary conditions have the form 

1 , .2n .27T . , . 

Wn\X) = —= exp(« — nx + t—my), (94) 
Vab a b 



and 



( T n)2 + ( T m)2 (95) 



for all (positive and negative) integers n and m. 

As |^(^)| 2 = 1/ab for all levels Eq. ( |T2"D which determines energy levels after the intro- 
duction of a ^-function potential ([11]) takes the form 

v ^wh: = 1 (96) 

n " 

with v = X/ab. 

Unperturbed eigenvalues have multiplicity 4 (for non-zero m, n) due to the existence of 
positive and negative values of m, n. To remove this degeneracy we consider in the above 
sum only positive integers and to have the same mean density (p = ab/An) we divide all 
eigenvalues by 4 after which eigenvalues included in the sum are 

en = {-nf + Crmf (97) 
a b 

and m, n > 0. 

The sum (^6]) formally diverges and for computation we consider the following renormal- 
ization 

V ^ 1 pEmax ^ 1 ^ ( ) 

p ft E e n J E m i n E e 

where E min and E max are minimal and maximal values of energy included in the sum. The 
subtracted integral (considered as principal value) equals log(E max — E)/(E — E min ) and one 
obtains the same relation between bare and renormalized coupling constants as before (cf. 

^ = -+iog|Ny^. (99) 

v ' pv E - E min 

We take v 1 = 1 and compute 100000 energy levels for such model. In Fig. |1] the cumulative 
nearest-neighbor distribution of these levels, N(s), is presented. This quantity equals the 
integral over the nearest-neighbor distribution 
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Ms) = / p(s')ds' (100) 
Jo 

and it is better defined numerically than the usual nearest-neighbor distribution. On the 
same figure two other curves are presented. The dashed line corresponds to the Wigner sur- 
mise for the cumulative nearest-neighbor distribution in the Gaussian Orthogonal Ensemble 
(GOE) of random matrices 0| 

N GOE (s) = l-e-™ 2 l\ (101) 
The thin solid line represents the cumulative nearest-neighbor distribution for the so-called 



semi-Poisson model [12], |13] which serves as a reference point in models with intermediate 
statistics 

N s , P (s) = 1 - (2s + l)e~ 2s . (102) 

It is clearly seen that the cumulative nearest-neighbor distribution for the Seba billiard is 
quite far from GOE result and it is in between the semi-Poisson curve and the GOE one. 

The numerically computed 2-point correlation function for this model is plotted in Fig. |2|. 
The two curves in this figure correspond to theoretical predictions for small and large values 
of the argument given by Eqs. (|80|) and (^) respectively. 



V. SERIES EXPANSIONS 

The above expressions are quite cumbersome. Therefore it is of interest to represent 
them as power expansions. We start with function G(x,y) defined in ( |193| ). It is convenient 
to define 

G(x,y)=ig(t,3), (103) 

where 

t — ix, s — iy. (104) 
Using the standard formula for the Bessel function 

^ 2 m!(m + n)! 



one gets 



-t) m ^ (-*)» ^ t+s 



-e 



m=Q m - n=m n - 



co f j.\m m—1 ( \n 

1-E^-E M^ 8 . (106) 



m =0 m - nT0 n[ 



Expanding the exponent leads 
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00 (—f)m(_„\n 

g(t,s) = l- E , , (107) 

n m— n+/ — 1 

^Kn)^^-!) 1 E C ? n("l) fc , (108) 

Z=0 A.-0 

and are the binomial coefficients. Only terms for which the upper limits in these sums 
are non-negative are included in the summation. But 

E^(-l) fc = (-l) L ^_ 1; (109) 

k=0 

therefore 

n 

R(m, n) = £ ^C^" 1 = (-l) m+n ^C" +m -i- (HO) 

i=o 

Finally we get 

g(t,s) = l + E JiC-v (HI) 



Using Eq. ( |194| ) one can show that 

oo j.r 



1 oo j.m „n 

J( ai ,a 2 ) = s + t + 2 E -t-tC^^, (112) 



7T2C<J ^ -i 777. n 

n,m>l 



where, as before, i = —ia 2 /uj and s = ia\/uj. 

The expansion of the pre-exponent factor in fl68|) can be simplified by the following 
identity ( Pf, p.32) 



OO 



j2 (z) = V (-l) m (2m + 2n)! z 
nK ' ^ m\{m + 2n)\{{m + n)\] 2y 2 ) ' 1 } 



j 2 (O + jKO = i + E f n L iw M m (H4) 

[m!] 3 (m + 1)! 



One gets 



Changing a 2 — ► — «2 we can rewrite Eq. (|90D in the form 

roo ^ 2 

r 2 (fi) = 1 + / daida 2 P(ai, a 2 ) exp(-Q(n + —)ai - Q(n -)a 2 ) 

Jo v' v' 

exp(2z(ai + a 2 ) + 2niQQ(ai, a 2 )) + c.c, (H5) 

where 

P(a u a 2 ) = 1 + E f L XT TT (-«i^) m , (H6) 
[m!] 3 (m + 1)! 



and 



Q( ai ,« 2 )= e ^ +n -fd-^;™- 2 - (H7) 



m!n! ^+ m " 2 ' 

n,m>l 
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VI. LIMITING BEHAVIOR 



The above formulas give the exact expressions for the 2-point correlation function for 
the problem considered but they are quite cumbersome and suitable mostly for numerical 
calculations. The most interesting information which one can extract from them is the 
behavior of the 2-point correlation function at small and large Q. The purpose of this 
section is to discuss methods which permit to find these asymptotics without knowledge of 
the exact solution. 

It is clear that in order to find the behavior of the 2-point correlation function in the 
limit u} — > it is necessary to consider only the case when three initial levels (which we 
shall denote e\, e 2 , 63) are close each to others and all other levels are far from this triple. 
In other words only three terms in Eq. (fj7p are big. In such a case Eq. (fj7|) which should 
determine the positions of two nearest levels can be approximated as follows 

111 

+ - + - = 0. (118) 



E — e% E — e2 E — e% 
(Note the absence of the coupling constant.) The solution of this equation is 

ei + e 2 + e 3 , 1 



Ei, 2 = g ±-\]e{ + el + el- e x t 2 ~ eie 3 - e 2 e 3 . (119) 

This expression is translationally invariant, therefore one can assume that e\ + e 2 + e?, = 
and the difference between two adjacent levels is 



AE=^ef + el-e 1 e 2 . (120) 

After the corresponding rescaling the 2-point correlation function at the limit Q — > takes 
the form 

r 2 (ft) = \j s ( n ~ A=\fe\ + 4~ eie 2 )cfe lC fe 2 . (121) 

The factor 1/2 comes from the restriction e\ < e 2 . Changing variables e\ = rcos#, e 2 = 
r sin 9 and performing the integral over r one gets 

r 8 (n) = |n / - — — - ft . (122) 

8 J 1 — sin o cos u 

The last integral equals 4n/ \/3 and finally we obtain that in the limit of small f2 

r a (fi) = ^-n, (123) 



which coincides with the result ( p9|) obtained above by different method. 

To compute the behavior of the 2-point correlation function at large Q it is convenient 
to use a method based on the usual trace formula 

Let us define 
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x 1 



= (124) 

where all are independent random variables as before. 

We need to calculate the density of levels Ej defined by the equation 

vG(Ej) = 1. (125) 

Formally this density can be expressed in the following way 

p(E) = -hm[G(E) + A i g(i - vG(E))\, (126) 

where the symbol Im[F(i?)] means the following limit 

lm[F(E)} = lim — (F(E + ie) - F(E - is)), (127) 

taken over positive e. 

The derivation of Eq.( 126 ) is simple. The function 1 — XG(E) has zeros at Ej and poles 
at efc, therefore 

^ mi =j:_i__ ? _i_. (128) 

The first term in ( |126| ) cancels the poles from unperturbed levels and the imaginary part 
produces 5-function singularity at the required positions. 
Let us denote 

G±(E) = G(E±ie), (129) 

where e is positive and e — ► 0. 

In calculating the mean values it is useful to take explicitly into account the mean values 
of G±(E). Using the relation 

1 P-TirrS(x), (130) 



x ± ie x 
one finds 

< G ± (E) >= + plog(^t|), (131) 

where p is the mean level density of non-perturbed states. Introducing 

g±(E)=G ± (E)-<G±(E)>, (132) 

one can write 

1 - vG ± (E) = (l-v< G ± (E) >)(1 - \g±(E)), (133) 
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where 



v' 



X = - n . • A . ( 134 ) 
p{± ± Z7TV ) 



and v' is the renormalized coupling constant the same as in Eq. (|3q) 



1 1 W - E 

± = ± + \ n Z fL, (135) 

v' pv W + E K J 



The density of states (ignoring small correction to the mean density of state as in Eq. (39) 
now will take the form 

P{E) = P ~ ^n i9+{E) ~ 9 ~ {E)) " i^ (1 ° g(1 " Xg+{E)) ~ l0g(1 " y 9-(m- (136) 

The 2-point correlation function is the mean value of the product of two such expressions at 
different energies. The computation of the mean value can be done in perturbation theory 
by expanding this expression into powers of g±(E) and using a formula 



<gl(E 1 )g m (E 2 )> « p J 



de 



(E 1 -e + ie) n (E 2 -e-ie) 

m+n-2 AAm+n-1 ^ + 



2m(-ir- 1 pC^ + \_ 2 —t-—(l + 0(u)), (137) 



where u = E\ — E 2 - Therefore one can organize the perturbation series in a series of inverse 
power of us. 

Taking into account the first terms in the expansion of the logarithm in the above ex- 
pression one gets 



p{E) = p - ^[(1 - X-^)9+(E) - (1 - \ m j^)9-{E)]. (138) 



At large ui 



2 



R 2 (u) =p 2 - —{\— g+ {E 1 )g-{E 2 ) + (E 1 ~ E 2 ) + c.c.) 

= P 2 + 2 ( 2^1/ m » ( 139 ) 

UJ z {7l z + 1/u") 

which agrees with Eq. (|92| ) derived from the general formula. 

We stress that the methods used in this Section are not restricted to particular cases 
considered. They also can be used in more general situations where exact solution are not 
available, e.g. for rank-two perturbations (2 short-range impurities) and similar problems. 



VII. GENERAL CASE 

In the previous Sections we considered the calculation of the 2-point correlation function 
under the assumption of equality of all residues. Here the generalization of these calculations 
to the case of different residues is presented. 
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When the residues are different one has instead of Eq. (^7[) the following equation 

1, (140) 



v r 
E — e- 

j=l & C 3 



and, consequently, instead of Eq. fl49|) one gets a more general relation 



R 2 (Ei, -E 2 ) 



00 da\da 2 
-00 4tt 2 



-j(ai+«2) 



A' 



XI r k9(rkaur k a 2 ) f(r j a u r j a 2 ) 



k=l j+k 
r k i r k2 i)i(r ki ai,r ki a 2 )'ll) 2 (r k2 a 1 ,r k2 a 2 ) I] /( r i a l ? r i«2), 



(141) 



kij^k 2 



where /(«i, a 2 ), g(ai, a 2 ), and ^(ai, a 2 ) are the same as in Eqs. (|50|). Repeating the same 
steps as in Section [IV] one gets exact expressions for the 2-point correlation function. The 
analog of Eq. ([76]) (which is convenient for calculation of the small-f2 series of the 2-point 
correlation function) has the following form 



da\da 2 
4tt 2 ' 



-2irQ<J(rai ,ra2)>-iSJw(ai +»2)/i>' 



X 



+i 



d 2 



< $(r«i, ra 2 )e 



ir(ai— 02) ^2 



9 



5 



< $(rai, ra 2 )e 



ir(«i— 02) 



da.1 da 2 
x < $(rai,ra 2 )re ir(Ql - a2) 



> 



> 



Instead of Eq. fl9D|) useful for large-fi asymptotics one obtains 

r 2 (fi) = 1 + { / da x / ^ 2e - 2 ^< J ( ra i' ra 2)>-^( Q i+ a 2)/« 

x[< rJ 2 (2r v / ^^)e 4r(Ql - a2) > 2 

+ < rJ 1 2 (2r v /-ttia 2 )e ir(ai - Q2) > 2 ] + c.c.} 

Here < f(r) > denotes the mean value over all residues 

<f(r) >=-^£/fa), 



(142) 



(143) 



(144) 



3=1 



functions J(ai, a 2 ) and $(o;i, a 2 ) are defined in Eqs. (|78|) and (|172|) , v, and t> ' are 'bare' and 
'renormalized' coupling constants (see fl3"T|) and (|38|) ) . 

As in Section [IV] it is of interest to compute the behavior of the 2-point correlation 
function at small and large energy difference. When Q — > the integration by parts as in 



Section ITV] leads to 



r 2 (fi) -> £IA, 



(145) 



where 
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A = [ da ^ 2 [l < $( TO i,ra 2 )r 2 e ir(ai - a2) >< $(ra 1; ra 2 )e ir{ai - a2) > 2 
J Att 2 

+ < $(rai, ra 2 )re ir{ai - a2) > 2 < $(roi, ra 2 )e ir(ai -° 2) > . (146) 
Using ( |172|) this triple sum is transformed to the form (cf. with Eq. (|83|)) 

71 N 1 f°° 
JV i,j,k=l Z J0 
x f°° ^L e i(n+r ]+ r k )(v+im/2 + ac ( 147 ) 

JO 77 

According to ( P-69Q the last integral equals mH^^r^ + rj + r fc )£) therefore 

27T 2 w 1 Z" 00 

^ = ^3 E + / ^JofaO^ofaO^ofaO.Mfa + r,- + r k )Z). (148) 

iV ij',Jfc=l 

Taking into account Eq. (|85l) and symmetrising the answer one gets 



7T 1 



(r< + r 3 - + r k f 



a = y 

Of course, when all residues are equal, A = 7rv3/2 as in (|89|) . 
When — ► 00 from Eq. (|143|) one gets 



(149) 



r 2 (fi) -> 1+ < r > 2 / doti / d a2 e- n<r> ^(l Ql l + l a2 D-^( ai+a2 ) + c.c. 

JO J-00 

2 

~ 1 + n^ + i/v'^y (150) 



which differs from Eq. fl92|) only by a suitable definition of the coupling constant. 

Note that Eq. ( |149| ) is valid only for non-zero values of the residues. Otherwise, the 
pre-factor A formally diverges. This divergence is a consequence of a simple fact that when 
some rj = there exist certain energy levels exactly equal unperturbed levels. Therefore 
the set of new energy levels consists of two parts. The first includes energy levels which 
are changed by the perturbation. Their correlation function is given by the formulas above 
where only non-zero residues are taken into account. The second part consists of energy 
levels which are not changed by the perturbation. Their correlation functions are the same 
as for the Poisson process and, in particular, they do not display level repulsion. As the 
cross-correlations between (a finite number) of old and new energy levels disappear when 
iV — > 00, the resulting statistics is a superposition of two independent distributions and, in 
general, it will not have level repulsion (i.e. i?2(e) 7^ when e — > 0). 

The above case is realized e.g. when a 5-function potential (|11~D is added to a rectangular 
billiard with the Dirichlet boundary conditions and the positions of the singular point (xq = 
(#0,2/0)) are commensurable with the corresponding sides of the rectangular (a and b). In 
this model unperturbed wave functions are determined by two integers, n and m, 
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2 7T 7T 

ipfi = —== sin(— nx) smi—my) (151) 
\/ab a b 



and the residues are 



If 



— sin 2 (— nxo) sin 2 ( — my ). (152) 
ab a b 



xo Pi Vo P2 

— =—, — = — 153 

a qi b q 2 

for co-prime integers pi and qi there exist only a finite number of different residues depending 
on values n mod q\ and m mod q 2 . In particular, when n is divisible by q± or m is divisible 
by q 2 , fa = 0. It means that for all these values of n and m wave functions and energy 
eigenvalues will not be changed by the perturbation and the resulting distribution (included 
all energy levels) will not describe level repulsion. 

Another interesting case corresponds to a model when all residues are also independent 
random variables with a probability d[x(r). If rj never take very small values (more precisely, 
the mean value of 1/ y/r is finite) the only modification is that mean value over residues, 
< f(r) >, should be taken over the given distribution i.e. instead of Eq. (|144j) one has to 
use 



< f(r) >= J f(r)dfi(r). (154) 
In particular the value of the pre-factor A is 



A =H d t i(r 1 )d f i(r 2 )d^r 3 )J {ri+r2 + r3) \ (155) 
6 J V rir 2 r 3 



But, if the probability of small values of residues is large, certain expansions should be 
modified. A natural example is e.g. the Seba billiard with the Dirichlet boundary conditions 
when ratios of the positions of the singularity to the corresponding sides (as in ( |153|) ) are 
non-commensurable irrational numbers. In this case r n defined in ( |152| ) are equivalent to 
random variables 

r * = ab Sin2 ^ Sin2 ( ^ 156 ' ) 

where angles 0i and (fi 2 are uniformly distributed between and n. 

Now the 2-point correlation function at small Q will differ from Eq. ( pL45| ) . Indeed, formal 
calculation of the pre-factor ( |155| ) shows that it diverges at small r and its leading behavior 
is the following 

A^-<r><-^> 2 . (157) 
2 



But for variable (15 



<r>= — — ( T sin 2 dcf)) 2 = \, (158) 
ir 2 ab Jo ab 
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and 



1 Vab. f n d(fi 2 Vab, 2 , . 

where 0o is a cut-off of the integration over 0. With logarithmic accuracy 0o is proportional 
to Q, <po —>■ Q/Qo, and when — ^ 

r 2 (n) ^ ^- 3 mn 4 (n/n ). (160) 

The results of numerical calculations of 100000 levels of the Seba billiard with Dirichlet 
boundary conditions (with irrational ratios ( [L53|) and v' = 1) are presented at Figs. |3| and 
f|. In Fig. |H the cumulative nearest-neighbor distribution is plotted. The dashed and thin 
solid lines are the same as in Fig. II], the dotted line corresponds to the Poisson distribution 

N P (s) = 1 - e~ s . (161) 

Note that the computed distribution is quite far from all standard examples. Though the 
resulting distribution is more close to the Poisson distribution than the one for the Seba 
billiard with periodic boundary conditions (see Fig. [I]) one can check that this difference 
will be present for all non-zero values of the coupling constant (and, in particular, when 
v' — > oo). 

The 2-point correlation function is shown in Fig. |j. The limiting behavior for small (see 
( P-60|) ) and large (see ( |150| )) values of arguments are also indicated for comparison by thick 
solid lines. The value of fl in Eq. (|160| ), Q = 52.25, has been obtained by fitting expression 
( |160| ) to numerical result for small Q. 

VIII. CONCLUSION 

We have computed analytically the 2-point correlation function for zeros of random 
meromorphic functions with large number of poles when these poles are independent random 
variables. It is demonstrated that the statistics of these zeros corresponds to a distribution 
with level repulsion which differs from known examples. The resulting distribution is not 
universal but depends on residues. 

A natural example, where energy levels obey such meromorphic equation, corresponds to 
a rectangular billiard with a short-range impurity and our results give the spectral statistics 
of these models. It is of interest that different boundary conditions give very different 
results. Even the asymptotic behavior for small energy difference is different (cf. Eqs. ( |30"1) 
and (|I60D). 

We also proposed methods which permit to find the behavior of the 2-point correlation 
function at small and large arguments without the knowledge of the exact solution. The 
methods can be applied for cases where exact solutions are not available. 

IX. APPENDIX. 

The purpose of this Appendix is to compute the main integral (|52]) 
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f W Ot\ 0L 2 

I(a 1 ,a 2 )= (l-exp(i— ))de. (162) 

J-w ti\ — e E 2 — e 

Let us first derive a few useful relations. 

' 9 9 \u \ ■ f°° c ai 1 • "2 s de 

-)J{ai,a 2 ) = too / exp(z— M— J 7-=; rj-= r, (lodj 



9a 1 <9a2 ' Ei — e E 2 — e (Ei — e)(E 2 — e) 

where 

w = Ei- E 2 . (164) 
Perform in this integral the change of variable 



Now 



£ '- c = TT? < 165 > 



E 2 -e = - (#i - e)(^ 2 - e) = - ^^ 2 , de = ^ rft. (166) 



,9 9 . r ai-a 2 , f°° ( .a 1 a 2 dt 

tt— )J(ai,a 2 ) = -«exp(« )/ exp(z— (i -)) — . (167) 

aai aa2 J-00 00 a\t t 



Hence 



The following standard integrals will be useful for us (see jl4l). When Im/i > and lm/3 2 fi < 


;/ B 2 

t v ~ l exp(i^(t - ^—))dt = 2p u e iun/2 K_ l ,(pfx), (168) 
2 t 

and when Im/x > and Im/3 2 /i > 

[°° t 1 " 1 exp(i^(t + —))dt = iir(3 u e~ i7Tu/2 H w u ({3n). (169) 
J 2 t 

Here K u (x) and H^\x) = J±(x) + iY±(x) are respectively the Macdonald and Hankel func- 
tions. 

Note that K (x) is a real function, therefore 

f°° f a i u ot 2 dt f°° ,.<*!(. az 2 ,sdt 

/ exp(z— (t -))—=/ exp z— [t -))— -c.c. 

j-00 ait t Jo 00 ait t 

J 0, if aia2 > 

I 2mJ (— ^ v /aia 2 )sgn(ai), if aia 2 < ' 

and Jq (x) is the Bessel function of zero order. Finally 



(170) 



(tt ■7r-)J(ai,a 2 ) = exp(i— — — )$(ai,a 2 ), (171) 

aai aa2 lo 



where 
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$(ai,a 2 ) = 2TTj ((—y/-a 1 a 2 )Q(-a 1 a 2 )sgn(ai), (172) 
a; 



Q(x) = 1 if x > and Q(x) = if x < 0. Note that 

<9 2 ,1 



$(ai, a 2 ) = -o$(a!i, a 2 ). (173) 



Using the same method one can prove the following set of equalities 

d 2 i a± — a 2 



J(a 1 ,a 2 ) — exp(i )$(ai,a 2 ). (174) 

U) OJ 

J(ai, a 2 ) = exp(z— — — )-^-$(«i, a 2 ). (175) 
oo aoii 

ot\ — a 2 d 

J(ai,a 2 ) = — exp(i ) — — $(a!,a 2 ). (176) 



<9 2 «i - a 2 <9 

— -J J(ai,a 2 ) = -exp(« ) — 

oa 2 aa 2 

<9 4 T . , 1/5 9 u ai — a 2 



Note that the differentiation of function $ produces the (^-functions coming from the factor 

9(-aia 2 )sgn(ai) = -(sgn(ai) - sgn(a 2 )). (178) 
From Eq. (fL75) it follows that the second derivative 

da\ 

is known. Therefore 



J(«i, a 2 ) = 0i(«i, a 2 ) (179) 



J(a 1 ,a 2 ) = J{0,a 2 ) H ^- h/ («i - y)<j>i(y, at 2 )dy, (180) 

aai jo 

where J(0, a 2 ) and dJ(0, a 2 )/dai are values of J and its first derivative at a\ = 0. 
In Section |TTT| it was demonstrated that 

J(0, a 2 ) = ir\a 2 \ 

and symmetrically 

J(a%, 0) = 7r|ai|. 

As the second derivatives equal zero when aia 2 > the expression of J in these regions 
which is continuous when crossing the «i and a 2 axis is 

J(ai, a 2 ) = 7r(ai + a 2 )sgn(a 2 ). (181) 

It is clear that the function J{a\ 1 a 2 ) is a continuous function but with discontinuous first 
derivatives. The values of these discontinuities follow from the above discontinuity of the 
function $. 
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Therefore in the region 

cci > 0, a 2 < 0, ai + a 2 < 

(which is the continuation of the lower left square ot\ < 0, a 2 < through the negative 0:2 
axis) the function J {0.1,0.2) should take the form 



roii 

J (0.1,0.2) = -tt(q;2 + oi) + 2iraiexp{-ia 2 /uj) + / {ai — y)(f) {y,a 2 )dy, (182) 

j 



where 



and 



io(ai, "2) = exp(i— — — ) 7 ^-0(a 1 , a 2 ), (183) 
uj 001 



2 

0(ai,Q;2) = 27rJ (— v / -«i«2) (184) 

coincides with the function $ but without the discontinuous factor. 

After integration by part and certain transformations one obtains that in all regions 

Ct 2 Q.\ 

J{a u a 2 ) = n{a 1 + a 2 )sgn(a;2) - 7r[i(ai + a 2 )G{ , — ) (185) 

UJ UJ 

+ («2^o(0 + W-^iol 2 Ji{^)) exp(? )](sgn(ai) - sgn(a 2 )), 

where £ = ^\f—a\a 2 and 

roo 

G{x,y) = e iy J {2y/yi)e lt dt. (186) 

Jx 

The function G(x, y) obeys the following relations 

dG{x,y) 



dx 



-J (Oe M,+tw , (187) 



dG{x,y) . fx 



dy \j y 

where £ = 2^/xy. To prove the second identity note that 

dG{x,y) 
dy 

But the integrand equals 



-Ji{C)e ix+iy , (188) 



iy ["{Uo + J'^r-y'dt. (189) 
Jx \ y 



-i^yJIfyJ^JW)^)), (190) 



and as Jq = —Ji one gets the above relation. 
The functional equation 



G{x, y) + G{y, x) = % + i J {£)e ix+iy (191) 
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can be proved by comparing the derivatives of both sides and noting that the integral ( |T4[ . 
p.50) 

poo 2 1 a 2 

/ J (at)e-^ tdt = — exp(-— ) (192) 
Jo (27) 47 

requires that G(0, y) = i. 

When y < x the above integral can be indefinitely taken by parts which leads to expansion 

00 

G(x,y) = ie i *+ i yJ2(- i v) n U0, (193) 

n=0 



where rj = yy/x and £ = 2y/xy. 

The above collection of formulae permits to find the expansion of G(x, y) for all values 
of its arguments. 

Another useful representation of our integral is 

J(ai, a 2 ) = vr(ai + a 2 )sgn(« 2 ) (194) 

d d 

- ™[i(y - x)G(x, y) + (x— - y—)G(x, y)](sgn(x) + sgn(y)), 
where x = —a2/uj,y = ai/ui. 
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FIGURES 






2 
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FIG. 2. The 2-point correlation function for the Seba billiard with the periodic boundary 
conditions. Solid lines correspond to the asymptotics (80) and (p2|) for small and large values of 
energy difference. 
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FIG. 3. Nearest-neighbor distribution for the Seba billiard with the Dirichlet boundary con- 
ditions. Dashed line is the GOE result. Thin line is the semi-Poisson curve. Dotted line is the 
Poisson prediction. 
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FIG. 4. The 2-point correlation function for the Seba billiard with the Dirichlet boundary 
conditions. Solid lines correspond to the asymptotics (160) and ( |150| ) for small and large values of 
energy difference. 
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